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Abstract. We address here the issue of congestion in the modeling of crowd 
motion, in the non-smooth framework: contacts between people are not an- 
ticipated and avoided, they actually occur, and they are explicitly taken into 
account in the model. We limit our approach to very basic principles in terms 
of behavior, to focus on the particular problems raised by the non-smooth cha- 
racter of the models. We consider that individuals tend to move according to a 
desired, or spontanous, velocity. We account for congestion by assuming that 
the evolution realizes at each time an instantaneous balance between indivi- 
dual tendencies and global constraints (overlapping is forbidden): the actual 
velocity is defined as the closest to the desired velocity among all admissible 
ones, in a least square sense. We develop those principles in the microscopic 
and macroscopic settings, and we present how the framework of Wasserstein 
distance between measures allows to recover the sweeping process nature of the 
problem on the macroscopic level, which makes it possible to obtain existence 
results in spite of the non-smooth character of the evolution process. Micro and 
macro approaches are compared, and we investigate the similarities together 
with deep differences of those two levels of description. 

1. Introduction. Congestion phenomena in population dynamics cover a wide 
range of mechanisms, which can be classified according to their stiffness : 

1. Soft congestion: as the distance between individuals becomes smaller, the 
behavior of a single person is affected by the presence of others; 

2. Hard congestion: actual contacts between individuals occur, and the overall 
motion is perturbed by the fact that two persons may not occupy the same 
place at the same time. 

Soft congestion models can be included in both microscopic and macroscopic 
models. In the first class there are discrete-space models like cellular automata- 
based models which constrain pedestrians to be located at squares of a fixed grid j7j 
HT1 Hr?I [551 IM] . We can also mention models based on networks as route choice 
models [51 IB] or queuing models |40[ W2\ . A lot of evacuation softwares rest on such 
models as for example buildingExodus in [32]. Some microscopic models are space 
continous as the social force model introduced in |35) and its forerunner which has 
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been proposed by [30]. In [35], people are identified to particles submitted to the 
laws of Newtonian mechanics. 

As for macroscopic models, soft congestion is also usually favored. Some of these 
models (see for example [18], [13] and [14]) own their origins in vehicular traffic 
models, and deal with the one-dimensional case. In higher dimension, similarly to 
the microscopic case, many models rely on social forces. One solution, used e.g. 
by Bellomo and Dogbe in [4] and [24], or Degond in [22], is to add to the equation 
satisfied by the velocity a repulsive term that forces people to avoid high density 
areas. Another possibility is to modify directly the velocity by adding a term that 
make people deviate from their preferred path as soon as they approach a crowded 
area. The modeling of the velocity is the key point of these strategies, see for 
example the work of Hughes [36] , [37] , Coscia [19] , or Piccoli [47] , [48] . 

We aim here at addressing the particular issues pertaining to hard congestion. 
In this spirit, we shall make very crude assumptions on the social aspects triggering 
crowd dynamics (see Section [6] for possible improvements on these aspects), and 
simply consider that, at any time, a spontaneous (or desired) velocity field is given, 
and that it is purely reptilian: everyone's priority is to achieve his own goals, without 
accounting for others. Actual behavior shall result from some sort of compromise 
between individual tendencies and congestion constraints. More precisely, we shall 
consider that the actual velocity is defined as the projection of the spontaneous 
one onto the set of feasible velocities (i.e. which do not lead to a violation of the 
non-overlapping constraint). 

Those basic principles can be applied to microscopic and macroscopic descriptions 
of the crowd. In the microscopic setting, the degrees of freedom are the positions 
of individuals (identified to rigid disks) , and the non-overlapping constraints can be 
written straightforwardly by prescribing a minimal value for the distance between 
centers. The problem takes the form of a differential inclusion which fits into the 
general framework of sweeping processes introduced by Moreau in the 70 's (see [44]). 
This approach has been extended more recently to non convex cases, namely the 
case of uniformly prox- regular sets in [16], in [55] and later in [17] ■ The perturbed 
(i.e. inhomogeneous) problem has been studied in [TO] EU El] [28] . As we shall see, 
it gives a natural framework for the microscopic model we propose. 

In the macroscopic setting, the crowd is seen as a population density and the 
non-overlapping constraint consists in prescribing a maximal value for this den- 
sity. Although both microscopic and macroscopic models express the same type 
of modeling assumptions, the mathematical structure of the macroscopic model is 
less obvious. Expressed in the Eulerian framework as a transport equation by the 
actual velocity field, which is defined as the projection of the spontaneous one onto 
the set of feasible velocities, it involves a conservative transport equation by a field 
whose regularity cannot be controlled a priori, so that classical results cannot be 
used. A first attempt to address those issues was proposed in [43], in the case 
where the spontaneous velocity is the gradient of a given dissatisfaction function 
(typically the distance to the exit for the evacuation of a building): the framework 
of optimal transportation allows to reformulate the problem as a gradient flow in 
the Wasserstein space of measures (see e.g. [61] for an overview of the theory of 
optimal transport, and [3] for the notion of gradient flow in this setting). We show 
here that the sweeping process framework, which does not rely on any assumption 
on the gradient nature of the forcing term, can be extended in the macroscopic 
situation. 
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Stemming from similar principles, both microscopic and macroscopic models 
present nevertheless deep differences. For instance we shall point out the fact that 
the notion of maximal density is not properly defined in the microscopic setting. 
Thus the macroscopic one cannot be expected to be obtained from the microscopic 
one by any homogenization process. 

We shall focus here on the particular issues related to hard congestion. In order 
to highlight the very problems it raises, we will favor a very basic form of the 
crowd motion model. In particular we shall disregard in the main part of this 
paper any social or strategical aspects of human behavior. We are aware of the 
very crude character of these assumptions, and we gathered in Section |6] other 
aspects of pedestrian behavior that could be included in our approach. Besides, 
we hope that the framework presented here, and the issues we raise in attempting 
to establish links between microscopic and macroscopic settings, may be fruitful in 
other contexts, e.g. micro-macro issues in granular flows or modeling of chemotactic 
motion of large populations of cells. 

2. Hard congestion models. 

2.1. Microscopic setting. We consider a population of N individuals, identified 
with rigid disks centered at q 1; q2, . . . , q^, with common radius r > (it can 
be extended straighfowardly to the so-called polydisperse case, i.e. with different 
radii). We denote by 

U(q) = (U 1 (q) 1 ...,U JV (q)) 

the generalized spontaneous velocity of the crowd (U^ is the velocity which the 
individual i would like to have in the absence of others). We assume here that the 
desired velocity of i depends on his location only, and that this dependence is the 
same for all individuals (see Section [6] for more general behavioral models): 

U i = U (q i ), 

where Uo is a given field. The set of feasible configurations is defined as 

K = {qe WL 2N , Aj(q) = |<fc - <b| - 2r > , j} ■ 

We disregard here obstacles (like walls, furniture) to alleviate notations, but they 
can be included in the definition of K straightforwardly. Non-overlapping is pre- 
served by prescribing that Dij may not decrease if it is 0, i.e. Z)^ = Gy(q) ■ q > 
where Gy(q) = V-Dy. It leads to define the set of feasible velocities as 

C q = {v e R 2N , Ai(q) = =► Gy(q) • v > } . (1) 
The model simply writes 

f - P C ,U(q) (2) 
where the projection is performed according to the euclidean norm over M. 2N . 

Saddle-point formulation. Problem (|2|) is self consistent and drives the evolu- 
tion process, as will be shown in Section [5] Yet, the saddle-point formulation of 
the projection problem sheds light on the underlying Darcy-like structure of the 
problem, and introduces a pressure field which will have a straightforward interpre- 
tation in terms of modeling, and which will have a continuous counterpart in the 
macroscopic setting. 
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At some given configuration q, let us denote by A the set of couples i < j, 

such that i and j are in contact. Constraints on the velocity can be written in a 
matrix form (we drop the dependence of matrix B upon the configuration) 

Bv < 0, 

where the rows of B correspond to active constraints : 

G y -v>0 (i.j)eA. 

The saddle-point formulation of the projection problem consists in finding (u,p) G 
R 2JV x (where N\ is the number of active constraints), such that 

r u + b* p = u 

(3) 

[ Bu < 0. 

supplemented by the complementarity condition 

p ■ Bu = 0. 

Pressure p^ can be seen as the interaction force between individuals i and j. 

2.2. Macroscopic setting. In the macroscopic setting the crowd is represented 
by a density p, which we shall assume to be supported within some domain (the 
room) f2. In order to alleviate the notations, we shall present the model in the case 
of a closed room (see the end of Section 13721 for some details on the way an exit door 
can be accounted for in this framework). If we set at 1 the saturation value, the set 
of feasible densities is 

K = j/? G L^R 2 ) , J p(x) dx = l, supp(p) C O, < p(x) < 1 for a.e. a-j . (4) 

The density is advected by the actual velocity u (macroscopic counterpart of the 
microscopic one dq/dt) 

dtp + v ■ ( P u) = o, 

and u is defined as 

u = ^c p U, 

where C p is the set of feasible velocities. It can be defined unformally as the set 
of all those velocities which have a non-negative divergence on the saturated zone 
[p = 1] (where the density may not increase). More precisely, it is defined by duality 

as 

C p = jv G L 2 (n) 2 7 J v • Vq < V 9 G Hl(Q) , q(x) = a.e. on [p < 1] j , (5) 

with Hl(Q) ={g£ff I (H), q > a.e. in fi}. 

Saddle-point formulation. The dual expression of C p induces a natural saddle- 
point formulation for the projection problem, based on a pressure field p: Find 
{u,p) G L 2 {fl) 2 x H^VL), where H x p {Vt) = {q e #1(0), 1 = a - e - on [P < !]}> such 
that 

= U 

(6) 

< for all q G H ^(Q), 
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with the complementarity condition: 

u • Vp ■ 



o 

3. Rigorous formalisms, well-posedness issues. 

3.1. Microscopic model : generalized sweeping processes. We address here 
the question of existence and uniqueness of solutions to 

^ = P c U(q). 

dt q y ' 

where t M> q(t) g IR 2 ^ describes the motion of the crowd, C q is the cone of feasible 
velocities defined by ([T]), and U(q) is the desired velocity. Before stating the main 
results of this section, let us describe the general framework into which it fits, namely 
the sweeping processes, introduced by Moreau in the late 70's (see [44]). The original 
setting was the following: consider V a Hilbert space and t 1 — > K(t) C V a path of 
convex closed sets in V, with some regularity in time (e.g. K(t) is continuous for the 
Hausdorff distance). He was interested in the evolution of a point q(t) subjected to 
remain in K{t). Assuming that the evolution process tends to minimize the norm 
of 5, one ends up with the following process 

f t e-N K{t) ( q (t)), 

where Nx(q) is the outward normal cone to K at q: 

N K {q) = {veV, 3a > , qe P K (q + av).} (7) 

This normal cone has been introduced in |15| and is called more precisely the proxi- 
mal normal cone. Note that, as K is convex, this normal cone identifies to dlx , the 
subdifferential of the indicatrix function of K(t). The following discrete process (so 
called catching up algorithm), which is used by Moreau to establish well-posedness 
of the problem, gives a clear idea of the evolution mechanism. Consider a time step 
t > 0, an initial condition q° € K(0), successive positions are built according to 

As K{t n ) is closed and convex, the projection is well-defined and contractant, which 
allows to establish convergence results for the sequence of piecewise constant solu- 
tions q T . It is clear that the recursive process extends straightforwardly to the case 
where the projection onto K is properly defined (i.e. single-valued) into its neigh- 
borhood. The finite-dimensional sets satisfying this property were introduced by 
Federer in (29] under the name of positively reached sets. Then, they were called 
p-convex sets by Canino in |12| and later proximally smooth sets by Clarke, Stern 
and Wolenski in [15] . The final name "uniformly prox-regular set" will be given by 
Rockafellar et al. in [4T)1 [50] (see Definition 13.11 below) . 

The crowd motion model differs slightly from the sweeping process, as the fea- 
sible set is fixed whereas the point q (which represents the whole crowd) tends to 
evolve according to some given velocity U. The same principles can be extended 
straightforwardly, in particular: assuming U is not too large, and K is uniformly 
prox-regular, then for a sufficiently small time step t, the catching up algorithm 
can be used to build discrete solutions. In the same manner, one may write the 
evolution process as 

^-%(q(t))3U(q(f)). (8) 
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Note that this new formulation is a direct consequence of ([2]): as Nk(<1) and C q 
are mutually polar, the identity operator can be written as / — -Pjvv(q) + ^"C"q 
(see [45]). Equivalence between both formulations is not obvious, as an identity 
has been replaced by an inclusion, yet we shall see that it holds true under some 
conditions. 

The catching up algorithm reads as follows 

Step 1 (prediction): q n+1 = q" +rU(q"); 
Step 2 (correction): q' 1+1 = P K (q' l+1 ) . 

Let us now express in a rigorous manner that this is indeed an algorithm for r 
sufficiently small, and that this strategy allows to build a solution to our problem. 
Let us first give a proper definition of prox- regularity (see Theorem 1.3 in |50|). 

Definition 3.1. A closed set K is said to be 77-prox-regular at q if there exists a 
neighborhood O of q such that for all q G dK fl O and v G Nx(q), with \v\ = 1, we 
have 

B(q + i]v,?]) f)K = 0, 

where Nx(q) is the proximal normal cone to K at q (defined by ([7])). It is said to 
be uniformly prox-regular with a constant 77 if it is r\ prox-regular at every point of 
its boundary. 

As for our microscopic model, it holds 

Proposition 1. K is uniformly prox-regular with 

/ 12 

7? - r y iv(jv-i)(jv+i)" 

Proof. First it can be shown that ={q€ M. 2N , > 0} is uniformly prox- 
regular. In this case, tools of differential geometry can be used to compute the 
constant of prox-regularity which corresponds to the smallest radius of curvature 
(i.e. the largest eigenvalue of Weingarten operator). After calculation, we show 
that is uniformly prox-regular with constant ry/2. 

One may wonder whether the intersection of such sets (which is the case for 
K) is uniformly prox-regular with a constant depending only on the constants of 
prox-regularity of the smooth sets. From a general point of view, this is wrong 
as illustrated in Figure [TJ Indeed, we have plotted in solid line the boundary of 
a set S which is the intersection of two identical disks' complements. This set is 
uniformly prox-regular but its constant of prox-regularity (equal to the radius of 
the disk plotted in dashed line) tends to zero when the disks' centers move away 
from each other. In this situation, the scalar product between the normal vectors 
n! and n 2 (see Fig. [5]) tends to -1, which suggests that the prox-regularity constant 
also depends on the angle between normal vectors. 




FIGURE 1 . Vanishing of the constant of prox-regularity. 
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FIGURE 2. Evolution of the angle between the vectors ni and ri2. 



Consequently the proof of the uniform prox-regularity of K rests on a good 
estimate of the angles between gradients of active constraints. The existence of 
such a constant r\ relies on the positive linearly independence of gradients of active 
constraints. More precisely it can be shown that there exists 7 > 1 such that for all 



E Ay|Gy(q)| < 7 



E 

(iJ)eA(q) 



XijGijiq) 



(9) 

(i,j)eA(q) 

It can be deduced from this inequality that K is 77-prox-regular with r\ = 
(a detailed proof of this result can be found in [42l [58]). Furthermore, the given 
upper bound is obtained by considering a configuration q of aligned disks and by 
computing the constant of local prox-regularity at q (see Proposition 3.18 in [60] 
for more details). □ 

Note that the prox-regularity coefficient of K degenerates (i.e. goes to 0) as the 
size of the disks goes to zero, for a constant total mass. Furthermore the prox- 
regularity coefficient of K depends on TV with TV 3 / 2 for N large enough (see [52]). 

Now the following result can be established: 



Theorem 3.2. Let U be Lipschitz and bounded, for all T > and all qo £ 

there exists a unique absolutely continuous solution 1 1 — > q(i) to 

^+N K (q)3V(q) a.e. in [0, T] 
q(o) = qo- 

Moreover, this solution satisfies the differential equation 



^ = P Cq (U(q)) a.e. in [0,T] 
q(0) = qo- 



The well-posedness of the differential inclusion can be proved thanks to results in 
[2T1 128| . Then Proposition 3.3 in [5] claims that the solution satisfies the following 
differential equation 

^+P JVK(q) (U(q))=U(q) 

which is equivalent to ^ since the cones iVjf(q) and C q are mutually polar. 

3.2. Macroscopic model. We are concerned here with the question of existence of 
solution to the macroscopic model: given a desired velocity field x 1— > U(x) defined 
in f2, given an initial density p°, find p(x,t) such that 

9 4 p + V-(pu) = 0, (10) 
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where the actual velocity field u verifies 

u = P Cp U, (11) 

and C p is the set of feasible velocities defined by <j5j) . Eq. [10] is meant in a weak 
sense: 

/ pd tV + I /pu-V<^+ f <p(0,-)p° = V^C([0,r)xl 2 ). 
o Jn Jo Jn Jn 

The overall problem can be written as a non-local and nonlinear transport equa- 
tion 

d t p + V -(pu(p)) = 0, 

where the notation u(p) expresses a dependence of u upon the whole density field p, 
through the projection of U onto C p . The nonsmooth character of this projection 
and the fact that the regularity of u is not controlled (the projection is performed 
in the L 2 sense), rule out the possibility to use standard tools to define solutions to 
this equation. One may wonder whether the catching-up approach, which proved 
successfull in the microscopic setting, can be followed in the present context. The 
prediction step is straightforward: one transports p according to the desired velo- 
city field U during a time step r > 0. As for the correction (or projection) step, 
it appears immediately that the standard eulerian way to measure distances bet- 
ween densities (by estimating a given norm of their difference) is not suitable. In 
the microscopic setting, the catching up approach proved successful because the 
difference between two configurations identifies with a displacement field, and this 
is due to the very Lagrangian description of individuals. Recovering this feature 
calls for a new way to measure differences between densities, more respectful of the 
Lagrangian description of the motion, and this way is the Wasserstein setting. 



Optimal transportation and Wasserstein distance. We assume here that fl 
is convex. Let t : S7 — > SI be a measurable mapping, and p 6 CP(f2) a probability 
measure. We say that t pushes forward p onto v (written t#[A = v) if 

p(t-\A))= v {A) 

for any measurable set A. For given measures po and pi, we denote by II(/zo,/ii) 
the set of transport maps between po and p\. The quadratic Wasserstein distance 
W2(a*o, Mi) is defined by 

W 2 (po,Pi) 2 = inf / |t(x) - x\ 2 dp (x). 
ter%o,/w) Ja 

Notice that this is a sloppy definition that holds for atomless measures, which is 
anyway the case we are interested in; for general measures one should pass through 
the so-called transport plans, which we do not want to introduce here, and we 
address the interested reader to [3] EH ■ However, this quantity W2 can be proven 
to be a distance on the space CP(r2) and it makes the space of probability measures 
a geodesic space, i.e. every pair of points is linked by a curve (which is in this case 
a curve of measures), such that the length of this curve exactly equals the distance 
between the points. In case in the minimization defining W2 there is an optimal 
transport map t (which is the case if po is absolutely continuous with respect to the 
Lebesgue measure £ ), then this geodesic curve is known explicitly and it is given 
by Pt ■= ((1 - t)id + tt) # p . 
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In the spirit of the structure that we are imposing on T, we define the outward 
normal cone to a subset K C 7 (we actually apply the definition of the strong 
Frechet subdifferential in [3] to the indicatrix function of K) by 

v e N K {p) <^> [ v(x) ■ (t(x) - x)dp{x) < o (||t - id|| L 2 (p) ) 
Jn 

for all t such that t#p £ K . It makes it possible to express the problem like we did 
in the microscopic setting (Eq. (JSJ): 



and u transports p according to ([TO]) . 

We may now write the catching up algorithm (which we consider for the time 
being as a theoretical tool): given an initial density p° and a time step r > 0, build 
p , . . . , p n according to 



where the projection is performed in the Wasserstein sense and the set K is the 
constrained set introduced in (Q| . As for the microscopic model, we must check that 
both steps are well defined (possibly under some restriction on the time step r), 
and that the obtained discrete solutions converge to a limit which can be identified 
as a solution to problem (p~0|) (TTTT) . 

We shall assume here two properties of U. First we require that it tends to keep 
people within the room fi, i.e. that for r sufficiently small (id + rU) maps fi onto 
O. Obviously, in the case of an emergency evacuation, with an open exit door, we 
will of course alleviate this assumption by allowing U to cross the exit line (we will 
discuss later on how to "catch up'' the possible part of the mass that exits £1 because 
of this). Secondly, we need it regular enough (Lipschitz continuous is sufficient), so 
that, still for r sufficiently small (id + rU) preserves the absolute continuity of the 
measure, i.e. p « L d implies (id + tU)^ p <K L d . This assumption being made, 
the prediction step is well defined for any measurable field U. 

The key point is the projection step, which differs significantly from the mi- 
croscopic situation. It is a variational problem in the space of measures, which 
consists in minimizing the distance to a fixed measure among elements of K . As 
for existence, standard compactness arguments may be applied, but the question of 
uniqueness is more delicate. 

First of all, let us rule out the possibility to establish uniform prox-regularity in 
the spirit of Definition 13.11 Consider the one-dimensional situation. A Dirac mass 
at zero projects onto the characteristic function of (—1/2, 1/2). More generally, a 
combination of Dirac masses 



(we assume that the distance between supports of any two of them is always larger 
than half the total mass carried by the couple, so that they do not interact) projects 
onto the sum of characteristic function of the flattened Dirac masses: 



u - N K {p) 3 U, 



(12) 




(id + tU)^ p n transport (prediction), 
Pk (p n+1 ) projection (correction) 



(13) 





The Wasserstein cost can be computed straightforwardly, it is 
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On the other way around, consider a density p £ K (i.e. p(x) < 1 almost everywhere) 
and u> the largest open set such that p{x) = 1 a.e. in to. This open set w is a 
countable union of open intervals of lengths (oin)i<n<N (with possibly ./V = +oo). 
Among all densities which projects onto K at p, the farthest is the combination 
of Dirac masses with weights ct n , supported at the middles of intervals, and it is 
at distance (J2 a l/ 12 ) 1/2 ■ As 

a consequence, r\ prox-regularity can be expected at 
some points, but r\ is not bounded away from 0. Note that some measures are the 
projection of themselves only, even if they saturate the constraint at every point 
of their support: consider e.g. a dense open set of total measure 1 in (—1, 1), and 
define p as the characteristic function of its complement. As the saturated zone 
contains no interval, it cannot be the projection of a measure which violates the 
constraint. 

Note that this non prox-regularity can be seen as a natural consequence of the 
fact that the property in the microscopic setting degenerates as the granularity gets 
finer (i.e. when the radius goes to 0). 

On the other hand, the set K enjoys some kind of convexity property: it is 
geodesically convex (see (3j) : 

Definition 3.3. K £ is said to be geodesically convex if for any po, p\ £ K, 

the geodesic curve pt joining po and pi belongs to K for all t £ [0, 1]. 

This property might suggest that uniqueness of the projection can be obtained 
by standard arguments: in a so-called Aleksandrov non-positively curved (NPC) 
length space (see [T)), considering 2 minimizers and the measure halfway along the 
geodesic, convexity properties of the distance function p M> W 2 (p,p) 2 lead to a 
contradiction. 

Unfortunately, as soon as d > 2, CP(51) is not NPC and the previous distance 
function even turns out to have concavity properties. 

Yet, despite the previous assertions (non prox-regularity and the concavity of the 
squared distance), other kind of interpolation curves (called generalized geodesies) 
can be used to assert well-posedness of the projection problem, without any restric- 
tion on the distance from K . 

Definition 3.4. (Generalized geodesies, see [3J, p. 207) A generalised geodesic 
joining po to pi with base p is a curve defined by 

p t = ((1 -t)r +£r x ) #/ o, 

where ro (resp. ri) is an optimal transport map from p to p$ (resp. p\). 

Notice that when ro = id this gives again the expression of a geodesic in 3> 
according to the distance W 2 ■ 

This generalized notion can be used to establish the following proposition 

Proposition 2. For any p £ V^l), the Wasserstein distance to the set K of ad- 
missible densities defined in ([4]) is attained at a unique point PkP- The projection 
operator Pk is continuous. 

Proof. It is easy to check (but one can refer to Lemma 9.2.1 p. 206 in [3]), that 
the square Wasserstein distance enjoys strict convexity properties along generalized 
geodesies, i.e. W 2 (p 7 p t ) 2 < (1 — t)W 2 (p, po) 2 + tW 2 (p, pi) 2 . This gives the unique- 
ness once we know that the set K is still convex along these generalized geodesies, 
i.e. that po,pi £ K implies p t £ K. This is obtained by a computation of the 
Jacobian factor of the map (1 — i)ro + iri, see [53]. □ 
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Now that the catching-up algorithm is defined properly, the obtained sequence 
of discrete solutions allows to obtain an existence result: 

Theorem 3.5. Let p T be the piecewise constant interpolation of the discrete densi- 
ties given by the catching-up algorithm h!3\) . IfXJ is a C 1 velocity field, and p° £ K , 
then p T converges as r tends to to a solution of the macroscopic problem 

r d tP +v-( P M) = o, 

I u = P Cp U. 

Proof. Let us first describe the discrete quantities we obtain thanks to the catching 
up algorithm. We denote by r™ +1 the optimal transport between the projected 
density p n+1 = P K p n+1 and and we write t™ +1 = (id + rU)" 1 (see figured]). 

These transport maps allow to define a discrete velocity as follows: 

v „ +1 = id-t^or^ 

T 




FIGURE 3. Definition of the discrete transport maps 

We then define two different interpolations of these quantities. First, the piece- 
wise constant interpolation is given by : 

I if t€]»r,(n + l)r]. 

Using [43] for the projection part, it is possible to prove that v T satisfies the following 
discrete decomposition 

U = v r + Vp T + e T , 

where e T converges uniformly to as r tends to 0. We also define a continuous 
interpolation of (p n ) n as follows 

p T (t, .) = (T? +1 ) #/9 "+\ if t G ]nr, (n + l)r], 

where T" +1 = ((t — nr)v™ +1 +t™ +1 o r™ +1 ). This second interpolation satisfies 
the transport equation at velocity v n+1 = v n+1 o (T" + 

It is possible to prove a priori estimates on these interpolated curves, and there- 
fore prove that they both converge to the same limit. The continuous interpolation 
of (p n ) n gives that the limit density satisfies a transport equation, and the discrete 
decomposition of the piecewise constant interpolation proves that the limit velocity 
is indeed the projection of the desired velocity onto C p . □ 

We finish this section by considering the case where, due to modeling reasons, 
we allow the vector field U to let the mass exit through a part of the boundary. 
This means that id + rU is no longer supposed to map f2 into f2 but the segment 
connecting x to x + t\J(x) is allowed to cross d£l in a prescribed subset T C d£l 
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which stands for the exit. In such a case the transport step of the algorithm does 
not need to be changed, but we have to face, during the projection step, a density 
p n+1 which is no longer supported in ft. To bring back the mass to the original 
domain we want to consider a modified set K, that we call Ky'. 

K r = {pe 9(0.) : p = pn+pr, supp(p r ) C T , < p Q {x) < 1 for a.e. x} . (14) 

The reason for this choice is the following: we consider that as soon as a particle 
reaches T, instead of following its movement after T, we leave it on it. This is done 
for simplicity, but it only means that we are no longer concerned with what happens 
to the particles that have reached T, not that they are really blocked on the exit. 
Obviously, we needed to withdraw the density constraint on T, so as to let particles 
stay on it, and also to represent the fact that T stands actually for everything that 
happens at the door and beyond. 

Once the set Ky is introduced, the projection step is done by defining 



which means that we take a measure whose support may go beyond f2 and we 
project it onto the measures over Q satisfying some extra density constraint. We 
stress that the same kind of argument (projecting on a set of measures concentrated 
on J7) could also be used, in the case with no exit, to withdraw the (yet, natural) 
assumption that id + tU maps into fi. 

The mathematical problem with this modified Ky is much trickier. One of the 
difficulties, that prevent the usual theory to be applied, is the fact that the set 
Ky loses some of the properties that K had previously (in particular it is no more 
geodesically convex: the geodesic - for the W2 distance - between two points of 
Ky could go out of Ky). In particular we have no clue about the uniqueness of 
the projection, but the algorithm works in the same way if we accept to take any 
minimizer of the distance. 

3.3. Gradient flow setting. In case the spontaneous velocity is the gradient of 
some dissatisfaction function (e.g. distance to the exit in case of an emergency 
evacuation), the microscopic model can be put into a gradient flow form: 



The function ^ that we consider is the total dissatisfaction, i.e. &(q) = J2i D(qi), 
where D(x) may be, as we said, the distance to the exit. Indeed, as K is uniformly 
prox-regular, the proximal normal cone identifies with the Frechet subdifferential of 
Ik hence dp = + iV#-(q) (see [5TJ for more details). 

In the macroscopic case, under the same assumption that U is a gradient -VD. 
then one defines the global dissatisfaction function in a continuous setting: 



and the overall process can be seen as a gradient flow in the Wasserstein space, i.e. 
p is advected by u, 




— e -dtp, tp = (^ + I K ), 
where dtp is the Frechet subdifferential of tp defined by 



c^(q) = {v G R 2N , ^(q) - p(cO - v • (q - q) > o (|q - q|)} . 




d t p + V • (pu) = 0, 
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with 

u € —dp(p) for a.e. t, 
where p = ^ + Ik , and dp is defined in the following sense (see [5] ) : 

Definition 3.6. The strong Frechet subdifferential of a function tp at p is the set 
of fields u such that for all transport maps t, the following inequality holds true 

<fi(p)+ / u(x) ■ (t{x) - x)dp(x) < y>(t # p) + o(||t-id||£2( p )). 

Notice that the definition we gave before of Nk is nothing but a particular case of 
this one, when tp = Ik. 

The main advantage of this gradient-flow formulation is the fact that it allows 
to handle vector fields U less regular than what we need in the general case (where 
it is natural to require U to be Lipschitz continuous, i.e. D £ C 1 ' 1 ). For the whole 
theory on gradient flows, first in a finite-dimensional setting, then in Hilbert spaces, 
and finally in metric spaces and particularly in J'(fi), we refer again to [3|. 

In particular, the catching up algorithms made by the coupling of a prediction 
and a correction step, may be replaced in the case of a gradient structure by a 
single-step procedure, called proximal algorithm, where 

q n+l g argmin |y(g) + l<? ~^ }■ 

This algorithm has also been formalized in a metric setting, where it is known as 
minimizing movements (see |211 [5]): it has first been applied to the case of 7(0,) 
with the Wasserstein metric by Jordan, Kinderlehrer and Otto in |38j . 

The details about a gradient flow approach to the macroscopic framework of 
crowd motion are contained in |43j . where both the case with no exit (i.e. with 
the constraint p € K) and with exit (p € Kr) are dealt with, the second one being 
much trickier than the first. 



4. Numerical solution. 

4.1. Microscopic model. Moreau's catching up algorithm suggests a strategy to 
build approximate solutions to the evolution problem. Yet, projecting a configu- 
ration q onto K is not straightforward. The scheme we propose extends ideas 
introduced in |41] for granular flows. It consists in performing a catching up step 
with K replaced by some kind of inner local convex approximation. More precisely, 
considering that the configuration q™ at time t n is known, the predicted configura- 
tion is obtained by 

q ,i+1 = q" + rU(q' 1 ). (15) 
Then q™ +1 is obtained by projecting q™ +1 onto 

q n+1 =iVq n+ \ (16) 

where K a stands for 



X q = {q , D tJ (q) + Gy ■ (q - q) > V* + j] 
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FIGURE 4. Theoretical and numerical projections. 



In FigurelU we illustrate the set K C R 2N , intersection of sets Kij whose boundaries 
are plotted in solid line. The set K^n is delimited by the dashed line. The theoretical 
and numerical projections, respectively q n+1 := P^(q™ + TU(q™)) and q™ +1 are 
represented (for two examples of U(q™)). Indeed, as K is uniformly prox-regular, 
the projection onto K of q n + rC/(q' 1 ) is well-defined for t small enough. The 
replacement of K by the convex set Kqn is convenient because it allows us to use 
classical numerical methods to compute this projection. Indeed, this projection 
problem can be reformulated in a saddle-point form which can be solved by Uzawa 
algorithm : find (q n+1 , A) E R 2N x (R+) "'"" 1 ' satisfying 

f q n+1 = q" +1 +5>y Gy(q") 

J Vz < j, D tj (q") + Gij (q») • (q" +1 - q») > (17) 

Note that the Lagrange multipliers are not unique in general. Indeed there exist 
configurations q™ with strictly more than 2N active constraints so that the associ- 
ated gradients Gy (q") are linearly dependent. 

The matrix appearing in Uzawa algorithm is C = B t B where B is the matrix 
whose rows are vectors Gy(q). In [42], the authors quantify how the condition 
number of matrix C varies with the parameter rjq (setting a lower bound of the 
local prox-regularity of K at point q) when C is non-singular. As a consequence 
we expect that the Uzawa algorithm converges less quickly for configurations with 
low local prox-regularity. In numerical simulations, we noticed indeed that solving 
the saddle-point problem requires more iterations in case of a jam. 

Numerical analysis. In the proposed scheme (|16|) . the replacement of K by the 
convex set Kqn is computationally convenient but arises many difficulties in the 
numerical analysis. The convergence result is based on the fact that Kq is a good 
approximation of K near q. 
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Theorem 4.1. Let us denote by q T the continuous piecewise linear function satis- 
fying q T (t™) = q" defined by Mb}) . Let U be Lipschitz and bounded, for all T > 0, 
the sequence (q T ) r uniformly converges to q defined in Theorem \3.°A when t tends 
to 0. 

Proof. Here we just give a sketch of the proof. We refer the reader to [58] for 
the details. It can be easily proved that (q T )r is bounded and so a convergent 
subsequence can be extracted. Thanks to the uniqueness result in Theorem 13.21 it 
suffices to check that the limit function (which is absolutely continuous) satisfies 
the differential inclusion. Since q" +1 = Pk^ (q" + rU(q")), it comes for all q, 

(q n + rU(q n ) - q" +1 ) • (q - q" +1 ) < |q" + rU(q n ) - q" +1 | d K ^ n (q). 

By dividing by r, we obtain for all q 

(U(q n ) - u n ) • (q - q" +1 ) < |U(q") - u"| d XqB (q) < 2\\V\Ud K ^ (q) 

where u" = (q™ +1 — q™)/r represents the discrete actual velocity. We aim at passing 
to the limit in the previous inequality and the crucial term is the last one. Thus 
the convergence result rests on the local continuity of the map (q, qo) — > d^^ (q) 
in a neighborhood of the set {(q, qo), q = qo}- More precisely it can be checked 
that for all q £ K and all q <E B(q, r q ), 

H n— >+oo 

by using the saddle point form of these projections (fTT)) . Thanks to the reverse 
triangle inequality, the Lagrange multipliers (which are not unique in general) are 
bounded. Then compactness arguments allow us to obtain the convergence of the 
projections. We deduce that the limit function q satisfies the following differential 
inclusion: 

^ + 7V(^ q ,q)9U(q) a.e. in [0,T]. 

As for all q G Q, N(Kq,q) = N(K,q), the required result is proved. □ 
The convergence order of this scheme will be specified in a forthcoming paper 
[B]. This more accurate result is based on metric qualification conditions between 
sets = {q, Dy(q) + G,j • (q — q) > 0}. More precisely there exists a constant 
9 > such for all q £ K and all q £ -B(q, r q ), 

d*,(q)<0XXy(q). 
This result rests on the reverse triangle inequality 

4.2. Macroscopic model. The strategy we propose relies again on the catching 
up philosophy which already made it possible to establish an existence result. The 
density is first advected by the spontaneous velocity field (with no consideration of 
the congestion constraint), and then projected (in the Wasserstein sense) onto the 
set of feasible densities. The time discretization scheme writes as follows: Given a 
time step r > 0, an initial configuration p°, approximate configurations p , . . . , p n 
are built recursively according to 

p n+1 = (id + rU) # p n (prediction step), (18) 
p n+i = p K pn+i (correction step). (19) 
The first step consists in transporting p n according to the given transport map 
id + rU. We simply transport the center of each cell at velocity U during a time 
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t, and then distribute the mass of the transported cell on the different cells of the 
mesh it intersects, as illustrated in Fig. [S] The transported density writes 

p n+1 = — ^- Pij A ( transported cell (i,j) n cell (k,l) ) , 

i,j k,l 

where A{B) represents the area of the set B. 
The second step is less standard. The scheme we propose to approximate the 
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Figure 5 . Lagrangian transport of the density 



projection onto K with respect to the Wasserstein distance is based on the following 
consideration^]: Consider a domain to and a non-zero density fx supported by uJ. 
For e > 0, l u + e\i ^ K, and its projection onto K is identically 1 in ui. Let us 
denote by su the displacement field which corresponds to the optimal transport 
between l w + sfx and Pk(1u + sfj,). One has 

(id + eu) # (l u + en) = Prt(1 w + en). 

As the latter is identically 1 in w, one has 

1 + = 1 

det(I + eVu) 

so that, at the first order in e, 

V • eu = Efi. 

As id + eu is an optimal map, u is a gradient : u = — Vp, where p verifies the 
Poisson problem 

-Ap = n- 

Since we are considering a violation of the density constraint on a set u> and assuming 
that, after the projection, the density will be saturated on oj itself, we are only 
interested in finding the mass that will exit u> when the displacement is given by 
eu. This means, at a first order approximation, we only need to estimate the flux of 



1 To alleviate the presentation, we do not normalize measures to a unit total mass, but the 
considerations on optimal transportation and Wasserstein distance for probability measures which 
we presented at the beginning of Section 13.21 are straightforwardly extended to measures with 
arbitrary mass. 
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density through the boundary du>, i.e. u- n = —dp/dn. Now consider the stochastic 
interpretation of the Poisson equation 

-Ap = ft. 

Considering that \i is a probability measure, we consider a random variable X G lo 
which follows the law of probability with density \i. We now consider a Brownian 
motion stemming from X. Its path crosses du for the first time at Y. The random 
variable Y <G du> is known to follow the probability law with density —dp/dn on dui 
(see [H]). 

The idea is therefore to redistribute the exceeding mass of the saturated zone in 
the following way (see Fig. [6]): For each saturated cell, a random walk is started, 
which transports the exceedind mass m = (p — 1)|C|, where |C| is the measure of 
the cell. When this random walk encounters a non-saturated cell, it gets rid of as 
much mass as it can, and continues as long as the transported mass is not fully 
distributed. When all the saturated cells have been treated, the obtained density 
p n+1 is admissible. 
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Figure 6 . Stochastic projection of the transported density 

Let us notice that the gradient flow structure when U = — V-D would also allow 
for another time-discretized algorithm, namely the proximal (or minimizing move- 
ment, or Jordan-Kinderlehrer-Otto) one instead of the prediction-correction one. 
This would lead to a sequence of minimization problems involving the computation 
of the distance Wi. The formulation by means of the transport plans, that we did 
not develop here for the sake of readability, transforms these problems into linear 
programming problems, that could a priori be solved through a simplex algorithm. 
Yet, it turns out to be so slow that 2D problems could not be efficiently solved in 
this way. 

Numerical analysis. The numerical analysis of the transport part of the algorithm 
is standard, but the projection part is more delicate, and we are not able to provide 
a rigorous error estimate for the stochastic scheme we propose. Let us simply say 
here that, beyond the unformal justification of the approach given previously, the 
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asymptotic behavior of the so-called Diffusion Limited Aggregation (DLA) process 
(see |39|). and the link between the distance based on the energy norm between 
measures (see [31] or [55]) and the Wasserstein distance give some arguments to 
support the chosen strategy. 

5. From micro to macro. Both microscopic and macroscopic model express the 
same assumptions at their respective levels: a spontaneous velocity is given, and the 
system evolves according to a velocity which is the closest to the spontaneous one 
among all feasible velocities, in a least square sense. Yet, the macroscopic model 
was not built from the microscopic one by a rigorous homogenization process. We 
describe here some obstacles to such a process, to shed light on the deep differences 
between both settings, in spite of formal similarities. 

Maximal density. 

First of all, let us point out that the notion of maximal density is somewhat am- 
biguous as far as the microscopic model is concerned. Let us consider the case of 
identical radii (monodisperse situation). The maximal packing density for identical 
disks is known to be p ma x = 7r/2v3 sa 0.9069 . . . , and corresponds to the triangular 
lattice (see Fig. [7]). Yet the actual density of moving collections of rigid disks is 




Figure 7. Triangular lattice 

generally strictly less than this maximal value, which is only attained for this very 
particular configuration. As an example, in the evacuation situation represented in 
Fig. [H the mean density upstream the exit door ranges typically between 0.85 and 
0.87. The fact that the actual density is strictly less than the maximal one does 
not mean that the flow is unconstrained (as the macroscopic setting would suggest). 
Those considerations call for a clear identification of configurations which saturate 
the constraint. Such a notion is proposed for the three-dimensional situation in |57| . 
in a Nash equilibrium spirit: 
We say that a particle (or a set of contacting particles) is jammed if it cannot be 
translated while fixing the positions of all of the other particles in the system. 
Note that this property is defined as local jamming in [57]. It is tempting to consider 
as maximal in some sense any density corresponding to such configurations, for 
which there are no free disks, so that constraints are activated everywhere. The 
triangular lattice is clearly jammed, but so is the cartesian lattice (p = 7r/4 rj 0.79), 
and it is possible to build looser jammed configuration (see Fig. [SI with p = 7r\/3/8 « 
0.68). 
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FIGURE 8. Maximal densities in the microscopic setting. 




Figure 9. Loose triangular lattice 



Constraints on the velocity. Beyond this fuzzy definition of maximal packing, 
the set of feasible velocities in the microscopic setting is strongly dependent on the 
local structure and not only on the density, a feature which is lost in the macroscopic 
approach. Let us first consider the triangular lattice represented in Fig. [7] Hav- 
ing the number of disks going to infinity and the radius going to zero, the density 
(characteristic function of the solid phase) converges to the uniform density p max ■ 
Microscopic feasible velocities can be defined in the solid phase, and one may won- 
der what are the corresponding feasible velocities in the macroscopic limit. Such 
velocities are clearly constrained in 3 directions. Using self-evident notations eo, 
e^ and e 27r /3 to design the principal directions of the lattice, one obtains at the 
limit some sort of unidirectional expansion constraints in those directions, which 
can be written 

e • Vu • e > , with e = eo , or e 27r /3. 

It still allows for non-rigid motions at the limit: Fig. [TBI represents a feasible velocity 
field at the microscopic level, whose macroscopic counterpart is 
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Note that, although the microscopic velocity field is tangent to the boundary of K, 
the corresponding macroscopic field is strictly expansive (V • u > 0). 




FIGURE f 0. Example of feasible velocity 



Jams. 

The considerations which we presented above have a crucial consequence on the 
behavior of both approaches: the microscopic model has the ability to reproduce 
blocked jams (which are observed in practice), whereas the macroscopic one does 
not. The latter property is a consequence of the maximum principle. Let us consider 
the situation represented in Fig. II H with a saturated zone (p = 1) upstream the 
exit door. The desired velocity U is assumed to point to the exit door, so that 
V • U < 0. The actual velocity is U — Vp, where p solves a Poisson problem 
in the saturated zone, with homogeneous Neumann condition on the walls, and 
homogeneous Dirichlet boundary conditions at the exit and on the interface with 
the non-saturated zone. By virtue of the maximum principle, p is nonnegative over 
the saturated zone, so that the velocity correction —dp/dn is non-negative on the 
door: people exit quicker than they would if there were no congestion (they are 
always pushed out by other people behind them). As a consequence, if the desired 
velocity field tends to move people out of the room, the evacuation process will 
never stop before the room is empty. 




Saturated zone 



dp/dn = p = dp/dn = 



FIGURE 11. Macroscopic evacuation 
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The microscopic situation is quite different. The pressures are still nonnegative 
by definition of the saddle-point formulation ([3]), but they are likely to act in the 
"wrong" direction for the individuals which are the closest to the door, as soon as 
they form an arch (see Figure [T2]) . 




Figure 12. Microscopic evacuation (arches) 



Mathematical issues (convergence or non-convergence). 

The precise way for setting convergence questions when passing from the micro- 
scopic to the macroscopic model needs to associate a density to every microscopic 
configuration. An easy way to do that is to associate to every non-overlapping disks 
configuration the uniform measure (with unit density) on the union of those disks. 
In this way we obtain a density p obviously satisfying p < 1 but we also know that, 
when we let the size of the disks go to zero (with the number of disks increasing 
consequently to infinity) we cannot obtain any density p < 1 in the macroscopic 
limit. In particular in dimension 2 it is impossible to go beyond 0.91, since we know 
that the maximal density is realized by the triangular lattice. 

This suggests that it is better to rescale the approximating p by the maximal 
density, which is 1 in dimension one, 7r/2v / 3 in dimension two. . . 

Now the question is the following: take a sequence of initial data p Q Nl i.e. the 
densities associated (in the way we described above) to a sequence of microscopic 
configurations with N non overlapping disks of radius r. Let N — > oo and N ~ r~ d 
and consider the limit density p°. Can we say that, for fixed velocity field U, the 
constrained evolution stemming in the microscopic model from p° N converges to that 
obtained in the macroscopic one from p°l 

The answer can be expected to be positive if d = 1, but the situation is much 
more complicated for d > 2. Actually, the unidimensional case gives no ambiguity 
between jammed and maximal density configuration, which is not the case in higher 
dimension. 

For instance one can consider a sequence of jammed cartesian lattices in a two- 
dimensional domain Q (a square, for instance), with a concentrating vector field 
U (for instance \J(x) = —x). In this situation, the configuration is completely 
blocked, and the evolution gives, for any time t > 0, the constant density p° N . On 
the other hand, the limit as TV — > oo is a constant density, but it does not activate 
the constraints. This is due to the fact that the density realized by the cartesian 
grid is strictly smaller than that of the triangular lattice, which was taken as a 
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reference. Hence, in the evolution, p t would differ from p and move towards a 
more concentrated configuration. 

But this (jammed configuration which are not of maximal density) is not the only 
source of problems in the micro-macro limit. Actually, one can consider a similar 
example with a sequence of triangular grids. In such a case, both p° N and the limit 
p activate and saturate the constraints. This induces some constraints on the 
admissible velocities but, as we saw in the paragraph devoted to these constraints, 
they act differently in the microscopic and in the macroscopic models. In particular, 
the limits of the admissible velocities for the microscopic problem are those vector 
fields which are "unilaterally incompressible" in the three main direction of the 
lattice, which is not at all the same as imposing only a positive divergence. The 
actual limit constrained space at the macroscopic level necessitates obviously to 
account for the microscopic structure of the contact network. Beyond packing 
fraction, different types of order metrics have been introduced in |25| to give a 
better description of the local structure (see also [57] for a recent review on this 
matter), but the use of those concepts for evolution problem is still widely open. 

This shows some deep differences, for d > 2, between the limit of the microscopic 
case and the macroscopic one. In the case of a gradient vector field U we could 
express this issue, thanks to the gradient-flow formulation, in terms of the associ- 
ated dissatisfaction functionals. Some general results give conditions to translate 
a variation form of convergence on these functionals (called T— convergence) into a 
convergence of the associated flows, but we do not want to insight this question in 
details because of its complexity; it would go out of the scope of this paper and it 
is moreover matter of an ongoing work. 



6. Possible extensions. 



6.1. Strategies. We specify here how the model we presented can be improved in 
order to account for more elaborate individual behaviors than the purely reptilian 
one on which we built our approach. 



Microscopic setting. We aim here at integrating the fact that individuals are 
likely to elaborate complex strategies to escape a building. For example, in case of 
congestion, they may decelerate or try to avoid the jam, instead of keeping pushing 
inefficiently The velocity of a person becomes then dependent upon the position 
of people he can see in front of him. Such a strategy can be defined as follows: We 



define the set Nj (see Fig. 13(a)) containing persons who are near and visible to 
the individual i: 

N 4 = { j, |q, - qj| < 2r + £ prox , dj ■ ey > cos a} , 

where dj = Uo(qi)/|Uo(q;)| an d ey = (q^ — qO/lqj — qi|- The constant a is taken 
equal to the half angle of view ( i.e. a ~ 60°). Two choices are possible for the 
individual i if his neighbours (belonging to Nj) walk slower than him. 

First, he can decelerate instead of going through the crowd. In this case, his 
speed s™ (i.e. his velocity's norm) at time t n is made dependent upon his neigh- 
bours' behavior at time t n ~ 1 . More precisely, it is computed as a barycenter of his 
neighbours' speeds, weighted by their relative positions. 

Otherwise, he can also be in a hurry and prefer changing his way instead of 
slowing down. In that case, if there exists another clear way through the set Nj, 
he follows it while keeping his desired speed, or else he will go round the group Nj. 
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(a) Illustration of the set Nj. (b) Bypass strategy. 



Figure 13. Integrating strategies. 

Among all directions allowing him to do so, he chooses d™ ew the closest to the one 
he wanted (see Fig. |13(b)| . 

Numerical simulations with these strategies are presented in [59] . If every pedes- 
trian prefers decelerating, there is typically no contact between people and rarefac- 
tion waves are observed. On the contrary, the avoiding strategy leads to situation 
where the constraints are activated. 



Macroscopic social forces. Integrating strategies in the macroscopic setting is 
more delicate, as it is less natural to follow the motion of a single individual in 
the crowd. However, it is possible to model the influence of local density on the 
behavior of the crowd by changing the desired velocity. More precisely, it is natural 
to assume that when people arrive upstream a crowded area, they tend to decrease 
their desired velocity. In the spirit of [361 |2"5] , a two dimensional generalization of 
road traffic models can be integrated in the model: at time t, people located at 
x € n have a desired velocity XJ p (x) with the same direction as the initial desired 
velocity U(x) and a norm that decreases when p increases: 

U p (x) = a(p)V(x), (20) 

where a(p) goes monotonically from 1 to as p goes from to 1 (saturation value). 
The evolution of the crowd then obeys the same principles as before, namely: 

d t p + w-{pu) = o 

u = P C V P - 

It can be proved that the density transported by U p stays in K: the projection 
step in the evolution equation becomes useless. Yet, if a(l) > (people continue to 
push even at saturation density), the congestion constraint is likely to be activated, 
and the approach we propose in this paper can be adapted to this situation. Note 
that, from a modeling standpoint, it is natural to use downstream information to 
determine the desired velocity (people adapt their velocity according to what they 
see, i.e. to the density in front of them, which is downstream their desired direction), 
whereas, once the velocity is determined, the transport can be performed using an 
standard upwind scheme. 

6.2. Multi-component populations. We investigate here the possibility to ac- 
count for different types of individuals, who might have different, and possibly 
antagonist, strategies. In the microscopic setting, because of the Lagrangian cha- 
racter of the approach, this can be done without any change in the mathematical 
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structure: it is sufficient to alleviate the assumption that the desired velocity de- 
pends on the position only, but can be defined differently depending of individuals. 
In the macroscopic setting, this is less straightforward, but it can be formulated in 
the following way: consider two populations 1 and 2, with associated densities pi 
and pi and desired velocity fields Ui and U2, we consider that in the saturated 
zone [pi + P2 = 1], desired velocities are perturbed by a common correction velocity 
w which ensures that the constraint p\ + P2 < 1 remains satisfied, and we define 
this velocity as the one that minimizes the L 2 norm. The model can be written 

" 9tPi + V-(pi(ui+w)) = 
dtP2 + V • (p 2 (u2 + w)) = 

(22) 

w + Vp = 
V • O1U1 + P2V2 + w) > 0. 

where p vanishes in the non-saturated zone [p\ + p2 < 1] , and verifies the comple- 
mentarity condition 

J (piU 1 + / 9 2 U 2 )-Vp = 0. 

Note that, in case Ui and U2 are both gradients of dissatisfaction functions (one for 
each population), the system presents a gradient- flow structure in the Wasserstein 
setting (see for a similar problem in the context of cell migration). 

6.3. Nash equilibrium approach. We give here some hints on another possibility 

to define the actual velocity, in a Nash equilibrium spirit. Note that this approach 

is likely to change the mathematical nature of the evolution problem. In particular, 

the instantaneous velocity is no longer defined in a unique way, but as any (among 

infinitely many) equilibrium in the following sense: The actual generalized velocity 

u = (m, U2, . . . , ujv) is such that, for any i, u considered as a function of u^ only 

minimizes the distance to U among feasible velocities, i.e. 

1 2 
u G argmin - |v - Uj| 
vec^(u) ^ 

where Cq(u) is the set 

C£(u) = {vel 2 , (ui,...,Ui_i,v,Ui + i,...,Uiv) G C q } . 

As soon as two disks are in contact and push against each other, there are infinitely 
many solutions to this problem, as illustrated by the situation of Fig. 03] the desired 
velocity of person 1 is 1 (in the horizontal direction) , whereas person 2 tends to stay 
still. In the previous approach (£ 2 projection onto the cone of feasible velocity), the 
actual velocity is 1/2 for both. In the present approach, any diagonal couple (a, a), 
with a G [0, 1], realizes an instantaneous Nash equilibrium. 




Figure 14. Two disks, Nash approach. 
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In the macroscopic framework, the problem can be described as follows: For any 
ui C M 2 , one prescribes that 

U| u G argmin - ||v| u - U| u || i2 , . 
vec^(u) * 

with 

C£(u) = {v G L 2 (n) , v © u lM 2 XlJ G C p } 

As indicated previously, this approach changes the mathematical nature of the 
evolution problem, which calls for further developments. Let us simply say that 
the solutions we built in the previous sections are particular solutions to this new 
problem, among infinitely many others. 



7. Numerical experiments. 



Micro-macro comparison. The first problem of the comparison between these 
models lies in the initial configuration: given an initial configuration for the disks 
in the microscopic setting, it is difficult to choose the macroscopic density that 
best fits this configuration. As long as uniform densities are considered, we adopt 
the following method: we first estimate the mean density in the initial microscopic 
configuration, and then normalize it with the maximal density of disks we encounter 
throughout the microscopic evolution. In Fig. 1151 we present an example of such a 
computation. 




Figure 15. Calculation of the macroscopic density: the initial 
density 0.53 (on the left) is normalized by the maximal density 
encountered 0.81 (on the right) to obtain a macroscopic density of 
0.65. 



We then let evolve both systems and compare the configurations at equivalent 
time steps. The formation of the saturated zones is quite the same in microscopic 
and macroscopic settings (see Fig. 116)1. However, as pointed out in Section [SJ the 
behavior of the two models at the exit is very different, in particular the evacuation 
is faster in the macroscopic model (see Fig. I17p. 




Figure 16. Formation of the saturated zones in the microscopic 
model (left) and the macroscopic one (right). 





Figure 17. Comparison of the evacuation for both models. 



Pressure field. In the microscopic model, the pressure exerted between people 
appears naturally as a Lagrange multiplier, and is computed using Usawa algorithm. 
The macroscopic pressure field does not appear explicitly in the numerical scheme, 
but it can be recovered by estimating during the projection step a discrete equivalent 
of the odometer function (see Section POj) : on each cell, we define the pressure as 
the total mass emitted by this cell during the stochastic projection. In Fig. [TBI we 
represented the microscopic pressure and the isovalues of the macroscopic pressure 
field upstream an exit door. In both settings, the pressure field is maximal in the 
middle of the saturated zone, and decreases at the entrance and the exit of this 
zone. 



Realistic geometries. Let us finally present some examples which correspond 
to more realistic evacuations situations. The microscopic setting, as suggested in 
Section [SJ describes more properly situations where people leave through narrow 
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exits. Fig.[19]shows the evacuation of the first floor of a large building (Departement 
of Mathematics at Orsay) . 

The macroscopic setting, on the other hand, would best fit situations where many 
people have to evacuate a large domain. We present in Fig.[20]the evacuation of the 
Stade de France. Even if the benches of the stadium are quite narrow passages, the 
macroscopic model gives results that correspond to real evacuation configurations. 

Acknowledgements. We would like to thank E. Maurel-Segala and L. Thibault 
for fruitful suggestions and comments. 
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FIGURE 18. Pressure field in the microscopic and macroscopic models. 
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FIGURE 19. Evacuation of the first floor of the Maths building at Orsay. 
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Figure 20. Evacuation of the Stade de France. 



